TEAM benchmark 7

https://www.compumag.org/wp/wp-content/uploads/2018/06/problem7.pdf

Solve Eddy-current problem

\[\begin{align*} \nabla\times(\nu\nabla\times\mathbf{A}) + j\omega\sigma \mathbf{A} = \mathbf{J} \end{align*}\]

with

\[\begin{align*} \mathbf{A}\times\mathbf{n} = \mathbf{0} \end{align*}\]

on the far boundary. The magnetic permability \(\mu_0 = 4\pi \cdot 10^{-7}\) H/m of vaccuum, and the electric conductivity of the aluminium plate is \(\sigma_{Al} = 35.26\cdot10^6\) S/m.

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
import matplotlib.pyplot as plt
%matplotlib widget
import numpy as np
# from netgen.meshing import MeshingParameters

Generate Geometry

The aluminum plate:

wpplate = WorkPlane(Axes((0,0,0),Z,X))
f = wpplate.Rectangle(0.294,0.294).MoveTo(0.018,0.018).Rectangle(0.108,0.108).Reverse().Face()
plate = f.Extrude(0.019)

plate.edges.hpref=1    # mark edges for geometric refinement
plate.faces.hpref=1    # mark edges for geometric refinement

plate.faces.maxh=0.04
plate.solids.name="plate"
plate.faces.col=(0.3,0.3,0.3,1)
Draw (plate);

The coil:

wp = WorkPlane(Axes(p=(0.294,0.000,0.049), n=Z, h=Y)) # origin in 

coil2di = wp.MoveTo(0.100,0.100).RectangleC(0.100,0.100).Offset(0.025).Face()
coil2di.edges.name="coili"
coil2do = wp.MoveTo(0.100,0.100).RectangleC(0.100,0.100).Offset(0.050).Face()
coil2d = coil2do-coil2di

cutout = wp.MoveTo(0.100,0.100).RectangleC(0.300,0.100).Face() + \
         wp.MoveTo(0.100,0.100).RectangleC(0.100,0.300).Face()

coil2d = Glue( [coil2d*cutout, coil2d-cutout])

# Draw (coil2d)
coil = coil2d.Extrude(0.100)
coil.solids.name="coil"
Draw (coil);

airbox = Box(Pnt(-.2, -.2, -.2), Pnt(0.5, 0.5, 0.5))
airbox.faces.name="outer"
airbox.solids.name="air"

airbox.faces.col = (0,0,1,0.3)

airbox = airbox-coil-plate
shape = Glue([ coil, plate, airbox])
Draw (shape);

geo = OCCGeometry(shape)

xi_msm = np.array([0, 18, 36, 54, 72, 90, 108, 126, 144, 162, 180, 198, 216, 234, 252, 270, 288]) * 1e-3
xi_sim = np.linspace(0,0.288,500)

mesh = Mesh(geo.GenerateMesh(maxh=0.2))
mesh.Curve(5)
mesh.RefineHP(levels=1, factor=0.2)
print("mats", set(mesh.GetMaterials()))
print("number of elements:", mesh.ne)

clipping_settings={"Clipping":{"enable":True, "z":-1, "dist":0}}
Draw(mesh, settings = clipping_settings);

mats {'air', 'plate', 'coil'}
number of elements: 19090

Coil currents

coilrect_xmin, coilrect_xmax = 0.294 - 0.150, 0.294-0.050
coilrect_ymin, coilrect_ymax = 0.050, 0.150

def Project(val, minval, maxval):
    return IfPos(val-minval, IfPos(val-maxval, maxval, val), minval)

projx = Project(x, coilrect_xmin, coilrect_xmax)
projy = Project(y, coilrect_ymin, coilrect_ymax)
tau_coil = CF( (projy-y, x-projx, 0) )
tau_coil /= Norm(tau_coil)
pot_coil = CF( (0, 0, sqrt((projy-y)**2+(projx-x)**2)-0.05) )

tau_coil_only = mesh.MaterialCF({"coil.*":tau_coil}, default=CF((0, 0, 0)))

clipping_settings={"Clipping":{"enable":True, "z":-1, "dist":-0.05}, "Objects":{"Clipping Plane":False, "Vectors":True}, "Vectors":{"grid_size":100}}
clipping_settings.update({"camera":{"transformations":[{"type":"move", "dir":(0,0,1), "dist":1.9}]}})
Draw(tau_coil_only, mesh,settings=clipping_settings, draw_surf=False); # top right, mesh, settings=clipping_settings)

FEM

mu = 4*pi*1e-7
sigma = mesh.MaterialCF({"plate":35.26e6}, default=1)

turns = 2742

def Solve(f):
    omega = 2*pi*f
    fes = HCurl(mesh, order=3, complex=True, dirichlet="outer", gradientdomains="plate")
    print("free dofs", sum(fes.FreeDofs()))
    u, v = fes.TnT()
    
    a = BilinearForm(fes, symmetric=True, condense=True)
    a += 1/mu*curl(u)*curl(v)*dx
    a += 1j*omega*sigma * u*v*dx

    # volume current density
    # f = LinearForm(-turns/(0.025 * 0.100)*tau_coil*v * dx("coil.*", bonus_intorder=4))

    # boundary current density + div-free correction
    f = LinearForm(
            -turns/(0.100)*tau_coil*v.Trace() * ds("coili.*", bonus_intorder=4) \
            +turns/(0.025 * 0.100)*pot_coil*curl(v) * dx("coil.*", bonus_intorder=4))
    
    A = GridFunction(fes)
    pre = Preconditioner(a, type="bddc", inverse="sparsecholesky")
    # a.Assemble()
    # f.Assemble()
    # inv = solvers.CGSolver(mat=a.mat, pre=pre, printrates='\r', maxiter=200)
    # A.vec[:] = inv*f.vec
    solvers.BVP(bf=a, lf=f, gf=A, pre=pre, tol=1e-12, solver_flags={"printrates": "\r"})
    
    B = curl(A)
    J = -1j*omega*sigma * A
    return {"A":A, "B":B, "J":J}

with TaskManager():
    ret = Solve(f=200)
B, J = ret["B"], ret["J"]

free dofs 460691

CG iteration 1, residual = 1.8783429375726386     
CG iteration 2, residual = 0.8770371762415737     
CG iteration 3, residual = 0.8567761852505471     

CG iteration 4, residual = 0.8672534502661857     
CG iteration 5, residual = 0.6706888994527486     
CG iteration 6, residual = 0.5418008998746271     

CG iteration 7, residual = 0.47862277985198426     
CG iteration 8, residual = 0.4922476312702463     
CG iteration 9, residual = 0.335523237543537     

CG iteration 10, residual = 0.40949399645394186     
CG iteration 11, residual = 0.2834753118364036     
CG iteration 12, residual = 0.31796900540176776     

CG iteration 13, residual = 0.24720542899765274     
CG iteration 14, residual = 0.24330692037309376     
CG iteration 15, residual = 0.047396792495629364     

CG iteration 16, residual = 0.05498441491280034     
CG iteration 17, residual = 0.1900991253071956     
CG iteration 18, residual = 0.09507346786932534     

CG iteration 19, residual = 0.02717188838096703     
CG iteration 20, residual = 0.023905530605082596     
CG iteration 21, residual = 0.044716229313266256     

CG iteration 22, residual = 0.4238120575843653     
CG iteration 23, residual = 0.06957774774464753     
CG iteration 24, residual = 0.03476951977409797     

CG iteration 25, residual = 0.02648196941574533     
CG iteration 26, residual = 0.04817581692516642     
CG iteration 27, residual = 0.019469444206568005     

CG iteration 28, residual = 0.011018584080742922     
CG iteration 29, residual = 0.008741504043485218     
CG iteration 30, residual = 0.006009965667606586     

CG iteration 31, residual = 0.008370870358348268     
CG iteration 32, residual = 0.008988519917875665     
CG iteration 33, residual = 0.016647146253052584     

CG iteration 34, residual = 0.022360496771420366     
CG iteration 35, residual = 0.008649210195311629     
CG iteration 36, residual = 0.005131236116501908     

CG iteration 37, residual = 0.002746893062508682     
CG iteration 38, residual = 0.005108037556395954     
CG iteration 39, residual = 0.004074821814958685     

CG iteration 40, residual = 0.007520869370844318     
CG iteration 41, residual = 0.0034470684291967434     
CG iteration 42, residual = 0.004331025626957288     

CG iteration 43, residual = 0.0023164800485316515     
CG iteration 44, residual = 0.01086682000882901     
CG iteration 45, residual = 0.0019416296733567923     

CG iteration 46, residual = 0.002237627940721444     
CG iteration 47, residual = 0.001366526055282683     
CG iteration 48, residual = 0.0020535076489096916     

CG iteration 49, residual = 0.001757128207996272     
CG iteration 50, residual = 0.0006804151469656269     
CG iteration 51, residual = 0.00032811138149604836     

CG iteration 52, residual = 0.00023387893442277193     
CG iteration 53, residual = 0.00036183717948650075     
CG iteration 54, residual = 0.0011531236803508736     

CG iteration 55, residual = 0.0006620220119144159     
CG iteration 56, residual = 0.0002282722285510288     
CG iteration 57, residual = 0.0005312076498532261     

CG iteration 58, residual = 0.001111217431539651     
CG iteration 59, residual = 0.0001309309701665514     
CG iteration 60, residual = 9.878574517357536e-05     

CG iteration 61, residual = 0.00031206590065498835     
CG iteration 62, residual = 0.00031368203472247335     
CG iteration 63, residual = 0.0002758615666564079     

CG iteration 64, residual = 0.0003290450223596543     
CG iteration 65, residual = 0.00019015340655358814     
CG iteration 66, residual = 0.0008149433857739929     

CG iteration 67, residual = 0.00017028842024570265     
CG iteration 68, residual = 0.0003443235333865015     
CG iteration 69, residual = 7.371359538943672e-05     

CG iteration 70, residual = 8.658172046253321e-05     
CG iteration 71, residual = 7.643645547007928e-05     
CG iteration 72, residual = 6.921864889655349e-05     

CG iteration 73, residual = 0.00011851056315627885     
CG iteration 74, residual = 0.00020579279776131612     
CG iteration 75, residual = 7.531176891107417e-05     

CG iteration 76, residual = 3.0483484436276334e-05     
CG iteration 77, residual = 2.1842848298331396e-05     
CG iteration 78, residual = 1.5860097509103003e-05     

CG iteration 79, residual = 1.920345679572137e-05     
CG iteration 80, residual = 3.0520691561658666e-05     
CG iteration 81, residual = 4.681928370927426e-05     

CG iteration 82, residual = 3.248975102233527e-05     
CG iteration 83, residual = 3.005671971446886e-05     
CG iteration 84, residual = 9.017468273393123e-06     

CG iteration 85, residual = 2.2164956667580047e-05     
CG iteration 86, residual = 1.0299063621415705e-05     
CG iteration 87, residual = 4.029542163018953e-05     

CG iteration 88, residual = 9.927399243983263e-06     
CG iteration 89, residual = 6.0378143114268436e-06     
CG iteration 90, residual = 4.20855490407563e-06     

CG iteration 91, residual = 1.3256804035461042e-05     
CG iteration 92, residual = 1.393260447538811e-05     
CG iteration 93, residual = 4.1256704649099325e-06     

CG iteration 94, residual = 8.547710465698398e-06     
CG iteration 95, residual = 2.198950696660927e-06     
CG iteration 96, residual = 1.2619219609978433e-06     

CG iteration 97, residual = 8.534512573574505e-06     
CG iteration 98, residual = 1.476733130397589e-06     
CG iteration 99, residual = 2.081384381260344e-06     

CG iteration 100, residual = 1.2847686130625406e-06     
CG iteration 101, residual = 8.765881777140543e-07     
CG iteration 102, residual = 3.0635170141866603e-06     

CG iteration 103, residual = 1.844956266482855e-06     
CG iteration 104, residual = 1.1512535791962352e-06     
CG iteration 105, residual = 2.292885224146807e-06     

CG iteration 106, residual = 9.828324876385927e-07     
CG iteration 107, residual = 8.2544265986808e-07     
CG iteration 108, residual = 1.55637974114929e-06     

CG iteration 109, residual = 5.510515618503845e-07     
CG iteration 110, residual = 2.65181428546799e-07     
CG iteration 111, residual = 4.854989945522712e-07     

CG iteration 112, residual = 8.485898438390512e-07     
CG iteration 113, residual = 8.143824244946461e-07     
CG iteration 114, residual = 8.667722184979637e-07     

CG iteration 115, residual = 7.356737568858367e-07     
CG iteration 116, residual = 1.2640638020527968e-07     
CG iteration 117, residual = 1.2463074640420159e-07     

CG iteration 118, residual = 3.458830122981088e-07     
CG iteration 119, residual = 3.4975623928864955e-07     
CG iteration 120, residual = 6.236916323613998e-08     

CG iteration 121, residual = 5.2143674565433173e-08     
CG iteration 122, residual = 2.8656701832876196e-07     
CG iteration 123, residual = 1.7062975345795272e-07     

CG iteration 124, residual = 7.494054920976881e-08     
CG iteration 125, residual = 1.326079341163922e-07     
CG iteration 126, residual = 3.956359788300151e-08     

CG iteration 127, residual = 3.495456133372768e-08     
CG iteration 128, residual = 1.2516222668111915e-08     
CG iteration 129, residual = 1.5831403798332296e-08     

CG iteration 130, residual = 2.481525826638752e-08     
CG iteration 131, residual = 2.2832552080249932e-08     
CG iteration 132, residual = 8.913970235940135e-08     

CG iteration 133, residual = 5.829724469955536e-08     
CG iteration 134, residual = 4.4172943334959026e-08     
CG iteration 135, residual = 3.061093396230331e-08     

CG iteration 136, residual = 2.4456324528578362e-08     
CG iteration 137, residual = 3.1020313166052286e-08     
CG iteration 138, residual = 1.5679565441881262e-08     

CG iteration 139, residual = 2.2336492509568013e-08     
CG iteration 140, residual = 7.826738435644444e-09     
CG iteration 141, residual = 1.0498845003452264e-08     

CG iteration 142, residual = 1.1057444383927013e-08     
CG iteration 143, residual = 4.592130187396944e-09     
CG iteration 144, residual = 6.319165765727571e-09     

CG iteration 145, residual = 3.94361782743871e-09     
CG iteration 146, residual = 1.0026433602850507e-08     
CG iteration 147, residual = 3.899198381688127e-09     

CG iteration 148, residual = 5.210956385031528e-09     
CG iteration 149, residual = 8.717933018740967e-09     
CG iteration 150, residual = 3.746093653942783e-08     

CG iteration 151, residual = 3.9000114650518835e-09     
CG iteration 152, residual = 4.530758238264773e-09     
CG iteration 153, residual = 5.6667314091816e-09     

CG iteration 154, residual = 3.14386634927285e-09     
CG iteration 155, residual = 2.374467476037341e-09     
CG iteration 156, residual = 9.343967173388642e-09     

CG iteration 157, residual = 4.036367475444364e-09     
CG iteration 158, residual = 6.697951493677635e-09     
CG iteration 159, residual = 3.096209850777817e-09     

CG iteration 160, residual = 1.5077771608774613e-09     
CG iteration 161, residual = 6.037473011281733e-10     
CG iteration 162, residual = 1.458456690004592e-09     

CG iteration 163, residual = 9.71027415822527e-10     
CG iteration 164, residual = 5.908798840475914e-10     
CG iteration 165, residual = 1.215922528829421e-09     

CG iteration 166, residual = 1.7195404722326009e-09     
CG iteration 167, residual = 6.88519621320981e-10     
CG iteration 168, residual = 1.2129413766462004e-09     

CG iteration 169, residual = 1.8214929303659361e-09     
CG iteration 170, residual = 3.23408467888459e-10     
CG iteration 171, residual = 6.410465326741992e-10     

CG iteration 172, residual = 3.7451756466943823e-10     
CG iteration 173, residual = 7.084887745234036e-10     
CG iteration 174, residual = 4.157296471179494e-10     

CG iteration 175, residual = 1.558672465415098e-10     
CG iteration 176, residual = 8.624885801209043e-11     
CG iteration 177, residual = 3.0515352535243603e-10     

CG iteration 178, residual = 1.3355917096631188e-10     
CG iteration 179, residual = 1.2074962717599965e-10     
CG iteration 180, residual = 1.0181083766898564e-09     

CG iteration 181, residual = 4.39574626846701e-10     
CG iteration 182, residual = 1.5989193067374444e-10     
CG iteration 183, residual = 1.9596446350420359e-10     

CG iteration 184, residual = 2.3112984341549125e-10     
CG iteration 185, residual = 3.3674317165139317e-10     
CG iteration 186, residual = 3.1536973010052365e-10     

CG iteration 187, residual = 3.052866562270005e-10     
CG iteration 188, residual = 1.5732874375987419e-10     
CG iteration 189, residual = 8.183288015159104e-11     

CG iteration 190, residual = 7.017746504985547e-11     
CG iteration 191, residual = 2.3702856682182534e-11     
CG iteration 192, residual = 8.219590526254736e-11     

CG iteration 193, residual = 1.6070540184121237e-09     
CG iteration 194, residual = 1.0607003319151e-10     
CG iteration 195, residual = 3.2165744513522936e-11     

CG iteration 196, residual = 1.5318806193615433e-11     
CG iteration 197, residual = 6.4901060854379464e-12     
CG iteration 198, residual = 1.1041055007166823e-11     

CG iteration 199, residual = 3.902368575630955e-11     
CG iteration 200, residual = 1.8076447671723485e-11     
WARNING: CG did not converge to TOL

Draw the B and J field

# clipping_settings={"Clipping":{"enable":True, "y":1, "z":0, "dist":0.1}, "Objects":{"Clipping Plane":False, "Vectors":True}, "Vectors":{"grid_size":100}}
# clipping_settings.update({"camera":{"transformations":[{"type":"rotateX", "angle":-90}, {"type":"move", "dir":(0,0,1.5), "dist":1}]}})
# Draw(B.real, mesh, settings=clipping_settings, max = 10e-3, min = 0, draw_surf=False);

clipping = { "function" : False,  "pnt" : (0,0.1,0), "vec" : (0,1,0) }
vectors = {"grid_size" : 50, "offset" : 0.5 }
Draw(B.real, mesh, clipping=clipping, vectors=vectors, max = 10e-3, min = 0, draw_surf=False);

# clipping_settings={"Clipping":{"enable":True, "y":0, "z":-1, "dist":-0.132}, "Objects":{"Clipping Plane":False, "Vectors":True}, "Vectors":{"grid_size":150}}
# clipping_settings.update({"camera":{"transformations":[{"type":"move", "dir":(0,0,1), "dist":2.2}]}})
# Draw(mesh.MaterialCF({"plate":J.imag}, default = CF((0, 0, 0))), mesh, settings = clipping_settings, min = 0, max = 8e5);

clipping = { "function" : False,  "pnt" : (0,0,0.0189), "vec" : (0,0,-1) }
vectors = {"grid_size" : 100 }
Draw(mesh.MaterialCF({"plate":J.imag}, default = CF((0, 0, 0))), 
     mesh, clipping=clipping, vectors=vectors, min = 0, max = 8e5);

Reference values of B and J from the bench-mark

Bz_A1B1_50Hz_ref = np.array([ -4.9  -1.16j, -17.88 +2.48j, -22.13 +4.15j, -20.19 +4.j,   -15.67 +3.07j,
   0.36 +2.31j,  43.64 +1.89j,  78.11 +4.97j,  71.55+12.61j , 60.44+14.15j,
  53.91+13.04j,  52.62+12.4j,  53.81+12.05j, 56.91+12.27j,  59.24+12.66j,
  52.78 +9.96j,  27.61 +2.26j])
               
Bz_A2B2_50Hz_ref = np.array([-1.83-1.63j, -8.5-0.6j, -13.6-0.43j, -15.21+0.11j, -14.48+1.26j, -5.62+3.4j,
 28.77+6.53j, 60.34+10.25j, 61.84+11.83j, 56.64+11.83j, 53.4+11.01j, 52.36+10.58j, 53.93+10.8j, 56.82+10.54j, 
 59.48+10.62j, 52.08+9.03j, 26.56+1.79j])

Jy_A3B3_50Hz_ref = np.array([0.249-0.629j,  0.685-0.873j,  0+0j,  0+0j,  0+0j,  0+0j,  0+0j,  -0.015-0.593j,
             -0.103-0.249j,  -0.061-0.101j,  -0.004-0.001j,  0.051+0.087j,  0.095+0.182j,  0.135+0.322j,  
             0.104+0.555j,  -0.321+0.822j,  -0.687+0.855j,])


Jy_A4B4_50Hz_ref = np.array([0.461-0.662j,  0.621-0.664j,  0+0j,  0+0j,  0+0j,  0+0j,  0+0j,  1.573-1.027j,
             0.556-0.757j,  0.237-0.364j,  0.097-0.149j,  -0.034+0.015j,  -0.157+0.154j,  -0.305+0.311j,
             -0.478+0.508j,  -0.66+0.747j,  -1.217+1.034j])

Bz on A1 to B1

Bz_A1B1_sim = np.array([B[2](mesh(x, 0.072, 0.034)) for x in xi_sim])

plt.figure(1)
plt.clf()
plt.title("Evaluate B.z on Line A1 - B1")
plt.plot(xi_sim, -1e4*Bz_A1B1_sim.real, "b", label="sim Bz.real")
plt.plot(xi_sim, 1e4*Bz_A1B1_sim.imag, "r", label="sim Bz.imag")
plt.plot(xi_msm, Bz_A1B1_50Hz_ref.real, "x", label="ref Bz.real")
plt.plot(xi_msm, Bz_A1B1_50Hz_ref.imag, "x", label="ref Bz.imag")
plt.ylabel("B in Gauss = $10^{-4}$ T")
plt.legend()
plt.show()

Bz on A2 to B2

Bz_A2B2_sim = np.array([B[2](mesh(x, 0.144, 0.034)) for x in xi_sim])

plt.figure(2)
plt.clf()
plt.title("Evaluate B.z on Line A2 - B2")
plt.plot(xi_sim, -1e4 * Bz_A2B2_sim.real, "b", label="sim Bz.real")
plt.plot(xi_sim, 1e4 * Bz_A2B2_sim.imag, "r", label="sim Bz.imag")
plt.plot(xi_msm, Bz_A2B2_50Hz_ref.real, "b--x", label="ref Bz.real")
plt.plot(xi_msm, Bz_A2B2_50Hz_ref.imag, "r--x", label="ref Bz.imag")
plt.ylabel("B in Gauss = $10^{-4}$ T")
plt.legend()
plt.show()

clipping = { "function" : True,  "pnt" : (0,0.072,0), "vec" : (0,1,0) }

Draw (J[1].real, mesh, clipping = clipping, min=-0.5e6, max=0.5e6 )
Draw (J[1].imag, mesh, clipping = clipping, min=-0.5e6, max=0.5e6 );

Jy_A3B3_sim = np.array([1e-6*1j*J[1](mesh(x, 0.072, 0.019-1e-5)) for x in xi_sim])
Jy_A4B4_sim = np.array([1e-6*1j*J[1](mesh(x, 0.072, 0.000+1e-5)) for x in xi_sim])

Jy on A3 to B3 - plate top:

plt.figure(3)
plt.clf()
plt.title("Evaluate J.y on Line A3 - B3")
plt.plot(xi_sim, Jy_A3B3_sim.real, "b", label="sim Jy.real")
plt.plot(xi_sim, Jy_A3B3_sim.imag, "r", label="sim Jy.imag")
plt.plot(xi_msm, Jy_A4B4_50Hz_ref.real, "bx", label="ref Jy.real")
plt.plot(xi_msm, Jy_A4B4_50Hz_ref.imag, "rx", label="ref Jy.imag")
plt.ylabel("J in $10^{6}$ A/mm$^2$")
plt.ylim([-1.5,1.5])
plt.legend()
plt.show()

Jy on A4 to B4 - plate bottom

plt.figure(4)
plt.clf()
plt.title("Evaluate J.y on Line A4 - B4")
plt.plot(xi_sim, Jy_A4B4_sim.real, "b", label="sim Jy.real")
plt.plot(xi_sim, Jy_A4B4_sim.imag, "r", label="sim Jy.imag")
plt.plot(xi_msm, Jy_A3B3_50Hz_ref.real, "bx", label="ref Jy.real")
plt.plot(xi_msm, Jy_A3B3_50Hz_ref.imag, "rx", label="ref Jy.imag")
plt.ylabel("J in $10^{6}$ A/mm$^2$")
plt.ylim([-1.2,1.2])
plt.legend()
plt.show()